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Let  X  be  a  positive  recurrent  regenerative  process  on  state  space  S 
with  steady-state  distribution  x.  Given  a  function  f  :  S  ♦  R,  we 
consider  the  problem  of  estimating  the  steady-state  central  moments 
p^(f )  “  /g  (f(x)-r)^  x(dx)  where  r  is  the  steady-state  mean  of 
f(X(«)).  We  obtain  strong  laws,  central  limit  theorems,  and  confidence 
Intervals  for  our  estimators,  and  present  numerical  results, 
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Let  X  ■  (X(t)  :  t  <  0}  be  a  (possibly)  delayed  regenerative  process 
taking  values  in  state  space  S,  with  regeneration  times  T(-l)  ■  0  <  T(0) 

<  T(l)  <  •••  (to  incorporate  regenerative  sequences  {XQ  :  n  .>  0),  we  pass 
to  the  continuous  time  process  X( • )  defined  by  X(t)  ■  X^j,  where  [t] 
is  the  greatest  Integer  less  than  or  equal  to  t).  Under  quite  general 
conditions  (see,  for  example,  p.  185  of  HEYMAN  and  SOBEL  (1982)),  there 
exists  a  probability  distribution  x  on  S  such  that 

(1.1)  r  (f)  m  4-  f  ds  ♦  /  f (y)  x(dy)  =  r(f) 

C  1  0  S 

a.s.  as  t  ♦  for  a  broad  class  of  functions  f:  S  ♦  R. 

An  important  problem  that  has  been  extensively  studied  in  the 

simulation  literature  concerns  the  estimation  of  the  parameter  r(f);  r(f) 

has  the  interpretation,  as  is  clear  from  (1.1),  as  the  steady-state  mean  of 

f(X(-)).  In  certain  applications,  however,  it  may  also  be  of  interest  to 

estimate  the  fluctuations  of  f(X(*))  around  its  steady-state  limit.  To 

be  precise,  set  *c(#)  "  f(#)  ~  r(f)  and  put  v(f)  ■  r(f^)  (for  g  :  S  ♦ 

define  g"  s  S  ♦  K  via  g"(x)  ■  g(x)  •  g(x)  •••  g(x)  (m  times)). 

2 

Letting  fc  play  the  role  of  f  in  (1.1),  we  observe  that  v(f)  may  be 
regarded  as  the  steady-state  variance  of  f(X(*)). 

More  generally,  let  u^(f )  •  r(f*);  then  p^(f)  is  the  m'th  steady- 
state  central  moment  of  f(*).  Clearly,  |ij(f).  ■  0  and  ”  v(f). 

Mote  that  estimation  of  u  (f)  "  r(fa)  is  not  a  special  case  of  the 

me 

standard  regenerative  method  (see,  for  example,  1GLEHART  (1978)),  since 
fc  depends  on  the  unknown  parameter  r(f),  which  Itself  must  be 


estimated.  Our  goal,  in  this  paper,  la  Co  develop  an  estimation 


■ethodology  for  the  central  moments,  In  Section  2,  we  prove  the 

required  limit  theorems  upon  which  our  methods  are  based*  Section  3 
develops  confidence  Intervals  for  central  moments,  and  numerical  results 
are  presented  in  Section  4. 


2.  ESTIMATORS  AID  LIMIT  THEOREMS  FOR  CENTRAL  MOMENTS 

He  assume  throughout  the  remainder  of  this  paper  that: 


(2.1)  ECYjCf2*)  +  xj“)  <  - 

where 


T(n) 

t  ■  T(n)  ~  T(n-l)  and  Y  (g)  ”  /  g(X(s))ds 
“  n  T(n-l) 


for  functions  g  :  S  ♦  JR.  Our  goal  is  to  estimate  the  central  moments 

p^f),  1  <  k  <m. 

The  following  binomial  representation  of  the  central  moments  is 
crucial  to  our  development: 


(2.2)  ^(f)  -  r(<f  -  r(f»)k 


-  I  (If)  r(fj)  (-I)11’*  r(f)k_J  . 
j-0  J 


For  g  :  S  •*  R,  set  r(n,g)  -  Y  (g)/x  ,  where 

n  n 


V*>  ■  i  J,  V«> 


and 


Relation  (2.2)  suggests  that 


(2.3) 


u(n,k)  -  l  (k)  r(n,  f^)  (-l)k“^  r(n,f) 

j-0  2 


k-j 


should  be  a  reasonable  estimator  for 


(2.4)  PROPOSITION.  u(n,k)  ♦  (i^(f)  a.s.  as  n  *►  »,  for  1  ^  k  2m. 


PBOOP.  From  (2.2)  and  (2.3),  it  is  clear  that  we  need  only  show  that 


k 

r(n,  f  )  ♦  r(f  )  a.s.  as  n  ♦  ■»,  for  k  <  2m.  The  strong  law  of  large 


m  lc 

numbers  applies  to  both  the  numerator  Y  (f  )  and  denominator  x  , 

n  n 


yielding  the  required  convergence,  provided  that  E|Y^(f  )|  <«  and 
ETj  <».  Clearly,  (2.1)  Implies  that  Etj  <  » ,  whereas  the  Inequality 


TO) 


E|Y1(fk)|  <  E(/  jf(X(s))|k  ds) 
1  T(0) 


T(l) 

<E(f  (1  +  f(X(s))ZB)  ds) 
T(0) 


.2«u 


-  Et}  +  EY^(f  )  <  • 


provides  the  flnlteness  of  the  other  moment. 


"•  *  - 
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By  Proposition  2.4,  u(n,k)  is  strongly  consistent  for  ^(f).  To 

~l/2 

state  our  next  result,  we  shall  use  the  notation  o^(n  )  to  denote 

1/2 

the  sequence  of  random  variables  (r.v. ’s)  6Q  such  that  n  6q  0  as 
n  ♦  *»,  where  *b  means  weak  convergence.  The  following  properties  follow 
easily  from  our  definition  and  standard  results  about  weak  convergence  (see 
p.  92  of  CHUNG  (1974)): 

(2.5)  i)  if  X  -  o  (n“1/2)  and  Y  ■  o  (n"l/2), 

n  p  n  p 

-1/2 

then  Z  -  X  +  Y  -  o(n  ) 
n  n  n  p 

-1/2 

ii)  if  1  4  1  and  T  •  o  (n  '  ), 
n  n  p 

-1/2 

then  Z  -  X  •  Y  -  o  (n  A,‘). 
n  n  n  p 

For  our  next  limit  result,  we  will  need  the  following  central  limit 
theorem  (CLT). 

(2.6)  ROPOSirSM.  For  l  <  k  <  m,  a1/2(r(n,fk)  -  r(fk))  -b  o(fk> 

N( 0 , 1 )  as  n  •  ,  where  o2(g)  -  EZ2(g)/(Et1)2,  Zn(g)  5  YQ(g)  -  r(g)TQ, 
and  N(0,1)  is  a  mean  zero  normal  r.v.  with  unit  variance. 

Proposition  2.6  is  well  known  in  the  regenerative  simulation 
literature  (see,  for  example,  Iglehart  (1978),  it  is  a  CLT  for  estimators 
of  the  uncentered  steady-state  moments  of  f(X(*)). 

(2.7)  PROPOSITION.  u(n,k)  -  r(n,f.)  +  o  (n"1/2)  for  1  <  k  <  m,  where 

K  p  —  — 

£k(*)  “  **<•>  -  kVl(£>  £c(,)* 
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PROOF.  Observe  chat 


(2.8)  u(n,k)  -  l  (k)  r(n,fJ)  (-l)k’J  r(f)k'j 

J-0  J 

+  V  (!f)  r(n,fJ)  (-l)k“J  (r(n,f)k_j  -  r(£)k"j)  . 

j-0  3 

Evidently,  for  0  £  p  £  k, 

(2.9)  n^2(r(n,f)*  -  r(f)*)  r(n,fp) 

-  nl^2(r(n,f)  -  r(f))  •  (  l  r(n,f)^  r(f)X  1  ^ )  r(h,fp)  . 

J-0 

From  the  proof  of  Proposition  2.4, 

(2.10)  (V  r(n,f)J  r(f)X"1-J)  •  r(n,fp)  *  JLr(f)l‘l  •  r(fp)  a.s. 

j-0 

as  n  -,  so  that  (2.9)  and  (2.10)  together  Imply  that 

(2.11)  (r(n,f)*  -  r(f)A)  •  r(n,fp) 

-  (r(n,f)  -  r(f))  •  lr(f)*  1  •  r(fp)  +  op(n  ^2)  . 

If 

By  noting  that  the  first  sum  on  the  right-hand  side  of  (2.8)  is  r(n,fc), 
we  can  combine  (2.8)  and  (2.11)  to  obtain 


u(n,k)  -  r(n,  f£)  +  (r(n,f)  -  r(f)) 

*  Y  (i)  r(fJ)  <-l)k"J  (k-J)  Kf)^"1  +  o  (n‘1/2) 

J-0  3  p 

-  r(n,  f*J)  -  r(n,  f  )  •  k 

c  c 

•  Y  C1"!1)  r(*j>  (-l)k"1_j  r(£)k‘J_1  +  o  <n"1/2) 

j-0  3  p 

-  r(n,  ffc)  +  op(n-1^2)  . 

Let  0(n)  -  (u(n,l),  u(n,m))T  and  |i  -  (pj(f) . |xm(f))T. 

He  take  all  vectors  as  column  vectors  and  T  denotes  transpose. 

(2.12)  THEOREM.  n1/2(0(n)  -  |»)  -»  N(0,  C(f))  where  N(0,  C(f)) 

is  a  multivariate  normal  r.v.  with  mean  vector  0,  and  covariance  matrix 

C(f)  with  elements  given  by 

C^Cf)  -  EZjCf^  Z^fjJ/CETj)2  . 

PROOF.  Let  a  be  a  arbitrary  column  vector  in  Ra,  and  note  that 
Proposition  2.7  Implies  that 


6 


>v- 


where 


n 

Z  (g)  s  n  l  Z.  (g)  ,  for  g  :  S  +  E  . 
n  lc-1  K 

In  the  second  equality  we  have  used  the  fact  that  “  rCf^).  Standard 

arguments  then  show  that 

n1/2(«T(U(n)  -  n))  -*  (aT  C(f)a)1/2  N(0,1) 

as  n  the  Cramer-Wold  device  (BILLINGSLEY  (1968),  p.  48)  then  yields 
the  desired  theorem. 

Theorem  2.12  shorn  that  our  estimators  have  an  asymptotically  norr  . 
distribution;  it  Is,  of  course,  a  limit  theorem  expressed  in  terms  of  an 
Index  n  corresponding  to  the  number  of  regenerative  cycles  simulated. 
However,  In  certain  settings.  It  Is  more  natural  to  express  limit  theorems 
in  terms  of  t,  the  amount  of  time  that  the  process  X  has  been 
simulated.  Then,  N(t)  ■  max(n  :  T(n)  <  t)  is  the  number  of  regenerative 
cycles  completed  by  time  t.  Set 


u(N(t),k);  N(t)  >  1 

u,.(k)  -  { 

z  0  ;  N(t)  <  1  , 

and 

Ot  -  (ut(l) . ut(m))  . 

(2.13)  PROPOSITION.  ut(k)  ^(f)  a.s.  as  t  -»•-»,  for  1  <  k  <  2m. 

This  follows  immediately  from  Proposition  2.4  and  the  fact  that 
N(t)  ♦  •  a.s.  as  t  ♦  *». 

(2.14)  THEOREM.  t1/2(0t  -  |i)  -4  N(0,  C*(f))  as  t  ♦  «,  where  C*(f) 
is  given  by 

C*j(f)  -  EZ1(f±)  Zj(f  )/Etj  . 

The  proof  of  Theorem  2.14  is  based  on  making  a  “random  time  change”  by 
substituting  the  process  N(t)  for  the  parameter  n  in  Theorem  2.12, 
giving 

N<t)1/2(U<N<t))  -  j»)  a  N(0,  C(f))  ; 

since  N(t)  a  t/Ex^,  we  obtain  the  result.  For  a  rigorous  proof  of  a  more 

general  result,  see  Theorem  3.10  of  GLYNN  and  IGLEHART  (1983). 

Specializing  our  results  to  the  steady-state  variance,  we  observe  that 

the  estimator  u(n,2)  (ut(2))  is  asymptotically  normal  with  limiting 

variance  given  by  EZ^^)2/^)2  (EZjCf^/EXj).  Since  EZ^)2  - 
2  2 

EZ.(f  )  ,  it  follows  that  the  asymptotic  variability  of  our  steady-state 


variance  estimator  is  unaffected  by  having  to  estimate  r(f).  (Note  that 
2  2  2 

EZ.(fc)  /(E tj)  is  the  variance  of  the  limiting  normal  r.v.  which 

2 

approximates  r(n,  fc>.)  For  higher-order  central  moments,  however,  the 
variances  will  generally  differ. 


3.  CONFIDENCE  INTERVAL  GENERATION 

In  this  section,  we  use  our  limit  theorems  of  Section  2  to  construct 
confidence  Intervals  for  the  k'th  central  moment.  To  accomplish  this,  we 
need  consistent  estimators  for  the  covariance  matrices  C(f)  and  C*(f). 


A.  (n.i.J)  -  /  (f(X(s))  -  r(n,f) )x  ds 

T(k-l) 


T(k) 

/  (f(X(u))  -  r(n,f))J  du  , 

T(k-l) 


(3.1)  PROPOSITION.  For  1  <  i,  J  <  m, 


An(i,J)  -  EY^f*)  YL(fJ) 


a.s.  as  n  » 


j<  V  >■  "a" 


*  -V-.-VV-. 

' '  ■_-  ‘  j  , 


(3.2)  PROPOSITION.  For  1  <  i,j  <  m,  C(n,i,j)  ♦  C^Cf)  a.s.  as  n  -*•  «. 

Application  of  the  converging-together  lemma  (BILLINGSLEY  (1968), 
p.  5)  to  Theorem  2.12  and  Proposition  3.2  shows  that  if  Ckk(f)  >  0 
(1  <  k  m) ,  then 

r  Ckk(f)1/2  CYkit)l/2, 

[u(n,k)  -  *(6)  - tjx —  ,  u(n,k)  +  z(6)  - jj^ - ] 

n  7  n  7 

is  an  asymptotic  100(1-6)2  confidence  interval  for  p^(f),  wtiere  z(6) 
solves  P(N(0,1)  <  z(6)}  *  1  -  6/2. 

A  similar  confidence  interval  can  be  based  on  simulation  of  X  to 
time  t.  Set 


CtUJ) 


C(N(t),  i.  j)  TN(t)  ;  N(t)  <  1 
0  ;  N( t)  <  1 


(3.*)  PROPOSITION.  For  1  <  i,  j  <  m,  Ct(i,j)  C*  (f)  a.s.  as  t  -*>  *». 


This  result  is  an  immediate  consequence  of  Proposition  3.2,  and  leads 
to  the  following  asymptotic  100(1-6)2  confidence  interval  for  p^(f) 
(assuming  (^(f)  >0.  1  <  k  <  m): 


c*k(f) 


1/2 


J72 


ut(k)  +  z(6) 


c*k<f) 

JT2 


1/2 


(3.5)  [ut(k)  -  s(6) 


For  the  steady-state  variance  v(f)  ■  ^(f),  we  can  use  our  knowledge 
that  iij(f)  ■  0  to  obtain  a  simpler  family  of  estimates  for  C 22 (*)  and 
CJ2(f).  Note  that 

(3.6)  EZj(f2)2  -  EZj(f 2)2 

-  EY.(f2)2  -  2v(f )  EY. (f  )x.  +  v(f)2  Ex?  . 
l  c  i  c  1  i 


C°(n)  -  ~2  {Aq(2,2)  - 


2u(n,2)  A  (2,0)  +  u(n,2)Z  A  (0,0)} 
n  n 


Cu(N(t))  x 


N(t)  ’ 


N(t)  >  1 


N(t)  <  1 


Propositions  2.4  and  3.1  show  that  <r(n)  ♦  C^j(f)  a.s.  as  n  ♦  ®  and 
C®  -*■  0^(0  a.s.  as  t  ♦  •,  provided  m  _>  2.  Thus,  if  C22(f)  >  0,  the 
following  Intervals  are  100(1-6)1  asymptotic  confidence  Intervals  for 
v(f): 


(3.7) 


r°f  \1/2  r°( 

[u(2,n)  -  z(6)  C-^2  ■,  u(2,n)  +  z(6)  — ] 
n  n 


(3.8) 


<C°)1/J 


«£>*« 


[utC2)  -  «(S)  -Sjtj-  ,  ut(2)  +  «<S)  t1y2  ] 


4.  NUMERICAL  RESULTS 


To  illustrate  the  results  obtained  in  the  previous  section  we  have 
simulated  three  models:  the  waiting  time  process  in  the  M/M/1  queue  (p  ■ 
0.5),  an  (s,S)  Inventory  model,  and  the  classical  repairman  model.  For  all 
three  models  we  selected  f  to  be  the  Identity  function  (f(x)  ■  x)  and 
k  •  2,  so  that  our  goal  was  to  estimate  the  variance  of  the  steady-state 
distribution. 

(4.1)  EXAMPLE.  M/M/1  Queue.  This  model  is  the  single  server  queue  with 

Poisson  arrivals  and  exponential  service  times.  We  simulate  the  waiting 

time  process  W  ■  (W  :  n  >  0} ,  where  W  is  the  waiting  time  (exclusive 

n  —  n 

of  service  time)  of  the  nth  customer.  Our  simulations  were  carried  out  for 

arrival  rate  X  -  5  and  service  rate  p  -  10,  and  hence  the  traffic 

intensity  p  -  0.5.  This  guarantees  that  Wq  W  as  n  >  «*. 

Regenerative  cycles  begin  at  those  values  of  n  for  which  Wq  •  0.  The 

2 

quantity  being  estimated  here  is  a  {W}  -  3.0.  We  did  50  replications  of 

5000  cycles  each.  The  sample  mean  of  the  50  point  estimates,  u(2,5000), 

1/2 

was  3.0417,  and  the  sample  mean  of  the  50  point  estimates  of  022(f)  , 

C®(5000)*^,  was  13.2782.  As  a  result,  the  sample  mean  of  the  90Z 
confidence  Intervals  was  [2.7328,  3.3506}.  The  coverage  probability  was 
58Z. 

(4.2)  EXAMPLE.  (a,S)  Inventory  Model.  This  model  is  a  periodic  review 
Inventory  model  with  a  stationary  (s,S)  ordering  policy.  An  (s,S)  policy 
is  characterised  by  two  positive  Integers:  s  and  S  with  s  <  S.  If  the 


•mount  of  Inventory  on  hand  plus  on  order  Is  lees  then  s,  order  to  bring 

the  sum  up  to  S.  If  the  inventory  is  greeter  then  or  equel  to  i,  do  not 

order.  Let  Xq  denote  the  level  of  Inventory  on  hend  plue  on  order  in 

period  n  immediately  after  the  ordering  dedalon.  If  d  denotes  the 

n 

demand  in  period  n,  then 


dn  —  *n 


s 


otherwise. 


We  assume  that  s  <  XQ  <  S.  The  state  space  of  {Xq  :  n  >  0}  is  {s,  s  + 
1 .  S>.  For  this  example  we  have  selected  s  ■  6,  S  ■  10,  and 

13/8  ,  J  -  0 

1/4  ,  j  -  1 

3/16  ,  J  -  2 

1/8  ,  j  -  3 

1/16  ,  j  -  4 

2 

Again  we  simulate  to  estimate  o  {X>  -  2.3333.  Using  i  -  10  as  the 
regenerative  state  we  run  SO  replications  of  1000  cycles  each.  The  sample 
mean  of  the  50  point  estimates,  u(2,1000),  was  2.33S2  and  the  sample  mean 
of  the  50  point  estimates  for  C22(f)1/2,  C°(1000)1/2,  was  1.2396.  The 
sample  mean  of  the  resulting  50  90%  confidence  intervals  was  [2.2708, 
2.3997].  The  coverage  probability  was  94X. 


(4.3)  EXAMPLE.  Classical  Repairman  Modal.  This  model  is  a  continuous 


time  Markov  chain  with  X(t)  denoting  the  number  of  failed  units  under¬ 
going  or  waiting  for  service  at  the  repair  facility  at  time  t.  We  have 


a  +  a  identical  aachines  each  with  an  exponential  failure  tiae  with 
failure  rate  X.  At  aost  n  of  the  unite  operate  at  one  tlae,  the  other 
a  being  thought  of  as  spares.  When  a  unit  falls,  it  is  sent  to  a  repair 
facility  consisting  of  s  repairmen  (servers)  having  exponential  repair 
(service)  tlaes  with  repair  rate  p.  With  these  assuaptlons  (X(t)  : 
t  £  0}  is  a  birth-death  process  with  state  space  (0,1,  ...  ,  a  +  n), 
birth  peraaeters  X^  •  (n  -  [i  -  a]+)X,  and  death  paraaeters  p^  •  p  • 
aln(l,s).  For  this  exaaple  we  have  used  n  -  10,  a  *  4,  X  ■  1,  p  ■  4,  and 
s  •  3.  Again  we  are  Interested  in  estiaating  the  steady  state  variance, 
o2{X)  ■  S.231.  We  ran  30  replications  of  1000  cycles  each  with  the 
regenerative  state  taken  to  be  1*2.  The  saaple  aeen  of  our  30  point 
estlaates,  u(2,1000),  was  3.1916,  and  the  saaple  aean  of  the  50  point 
estimates  of  C22(f)1/2,  C°(S000)1/2,  was  11.5562.  The  saaple  aean  of 
the  30  90Z  confidence  intervals  was  [4.5903,  5.7927].  The  coverage  prob¬ 
ability  was  74Z. 
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